# File combines data on energy use and GDP from WIOD 2016 release with exposures scores

library(readxl)
library(openxlsx)
library(readxlsb)
library(ggplot2)
library(concordance)
library(dplyr)
library(RColorBrewer)
library(gridExtra)

# Set working directory
path = "G:/My Drive/Research/AI and Energy/Drafting/ERL/Replication"
setwd(path)


################################################################################
################################################################################
# 1. Import and Prep WIOD Data
################################################################################
################################################################################
# Initialize dataframe to store relevant data
df.WIOD <- data.frame(matrix(nrow=0,ncol=5))
colnames(df.WIOD)=c("Year","Code","Description","Total Energy Use (TJ)","Output ($MM)")

# Get sector codes and description
df.categories <- as.data.frame(read_excel("Data/Input/WIOD/Gross10/NAMEA_GEU5610_USA.xls",sheet="Sector"))[c(1:57),]
df.categories$Code <- gsub("^.{0,1}", "", df.categories$Code) # Remove "v" in front of codes
# Get emissions data
df.emissions <- read.xlsx("Data/Input/WIOD/CO2 emissions.xlsx",sheet="USA")

# Loop over years and append yearly data to data.frame
for(yr in c(2000:2014)){
  print(yr)
  # Get energy data
  energy.yr <- read_excel("Data/Input/WIOD/Gross10/NAMEA_GEU5610_USA.xls",sheet=as.character(yr))[c(1:57),c(15)]$TOTAL
  # Get gdp data
  gdp.yr <- c(read_xlsb(paste0("Data/Input/WIOD/WIOT",as.character(yr),"_Nov16_ROW.xlsb"),range=paste0(as.character(yr),"!CYK2359:CYK2414"),col_names=F)$column.1,sum(read_xlsb(paste0("Data/Input/WIOD/WIOT",as.character(yr),"_Nov16_ROW.xlsb"),range=paste0(as.character(yr),"!CYK2359:CYK2414"),col_names=F)$column.1))
  # Get Emissions data
  emissions.yr <- c(df.emissions[c(1:56),as.character(yr)],sum(df.emissions[c(1:56),as.character(yr)]))
  df.tmp <- data.frame("Year"=yr,"Code"=df.categories[,1],"Description"=df.categories[,2],"Total Energy Use (TJ)"=energy.yr,"Output ($MM)"=gdp.yr,"CO2 Emissions (kt)"=emissions.yr)
  df.WIOD <- rbind(df.WIOD,df.tmp)
}
colnames(df.WIOD)=c("Year","Code","Description","Total Energy Use (TJ)","Output ($MM)","CO2 Emissions (kt)")


# Deflate prices to 2017USD
gdp.def = read.csv("Data/Input/A191RD3A086NBEA.csv")
gdp.def['Year'] = unlist(lapply(gdp.def['DATE'],FUN=function(x){as.numeric(substr(x,1,4))}))
gdp.def = gdp.def[,c(2,3)]
colnames(gdp.def) = c("GDP_Def","Year")
gdp.def['GDP_Def'] = gdp.def['GDP_Def']/100
df.WIOD = merge(df.WIOD,gdp.def,by="Year")
df.WIOD['Output ($MM)'] = df.WIOD['Output ($MM)']/df.WIOD['GDP_Def']

# Change units for analysis
df.WIOD['Output ($BB)'] = df.WIOD['Output ($MM)']*1e-3
df.WIOD['Total Energy Use (PJ)'] = df.WIOD['Total Energy Use (TJ)']*1e-3


# Create useful parameters
df.WIOD['Energy Intensity (PJ/$BB)'] = df.WIOD['Total Energy Use (PJ)']/df.WIOD['Output ($BB)']
df.WIOD['Emissions Intensity (kt/$BB)'] = df.WIOD['CO2 Emissions (kt)']/df.WIOD['Output ($BB)']
df.WIOD['Emissions Intensity (kt/PJ)'] = df.WIOD['CO2 Emissions (kt)']/df.WIOD['Total Energy Use (PJ)']




################################################################################
################################################################################
# 2. Make projected WIOD dataset going through 2023
################################################################################
################################################################################

# Log-linear projection output, energy, and emissions at industry level (Linear and quadratic projection leads to negative values)
yrs = c(2000:2023)
Code = unique(df.WIOD$Code)
Description = unique(df.WIOD$Description)
df.WIOD.proj = expand.grid("Year"=yrs,"Code"=Code)
df.WIOD.proj <- merge(df.WIOD.proj,df.WIOD[df.WIOD$Year==2014,c("Code","Description")],by="Code")
df.WIOD$Year2 = df.WIOD$Year^2
proj.cols <- colnames(df.WIOD)[4:6]
df.WIOD.proj[proj.cols] <- NA
for(i in proj.cols){
  for(j in Code){
    # lm_formula <- as.formula(paste0("`",i,"`"," ~ Year + Year2"))
    lm_formula <- as.formula(paste0("log(`",i,"`)"," ~ Year "))
    reg.data = df.WIOD[df.WIOD$Code==j,]
    reg.data[reg.data[,i]==0,i] <- 0.1
    regress <- summary(lm(lm_formula,data = reg.data))
    df.WIOD.proj[df.WIOD.proj$Code==j,i] = exp(regress$coefficients[1] + regress$coefficients[2]*yrs) #+ regress$coefficients[3]*yrs^2
  }
}

# Recalculate aggregate
for(i in proj.cols){
  for(y in yrs){
    df.WIOD.proj[df.WIOD.proj$Code=="TOTII_new" & df.WIOD.proj$Year==y,i] = sum(df.WIOD.proj[df.WIOD.proj$Code!="TOTII_new" & df.WIOD.proj$Year==y,i] )
  }
}


# Change units for analysis
df.WIOD.proj['Output ($BB)'] = df.WIOD.proj['Output ($MM)']*1e-3
df.WIOD.proj['Total Energy Use (PJ)'] = df.WIOD.proj['Total Energy Use (TJ)']*1e-3

# Create useful parameters
df.WIOD.proj['Energy Intensity (PJ/$BB)'] = df.WIOD.proj['Total Energy Use (PJ)']/df.WIOD.proj['Output ($BB)']
df.WIOD.proj['Emissions Intensity (kt/$BB)'] = df.WIOD.proj['CO2 Emissions (kt)']/df.WIOD.proj['Output ($BB)']
df.WIOD.proj['Emissions Intensity (kt/PJ)'] = df.WIOD.proj['CO2 Emissions (kt)']/df.WIOD.proj['Total Energy Use (PJ)']




pdf("TablesAndFigures/SI/WIOD_Projected/OutputProjection_Aggregate.pdf",width=6,height=5)
ggplot() + 
  geom_line(data=df.WIOD[df.WIOD$Code=="TOTII_new",],aes(Year,`Output ($BB)`,linetype="Data"),linewidth=1) +
  geom_line(data=df.WIOD.proj[df.WIOD.proj$Code=="TOTII_new",],aes(Year,`Output ($BB)`,linetype="Projection"),linewidth=1) +
  theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + xlab("Year") + ylab("Output ($BB)") + 
  ggtitle(paste0("GDP - ","Aggregate")) +
  scale_linetype_manual(name="Source",values=c("Data"="solid","Projection"="dashed")) 
dev.off()

pdf("TablesAndFigures/SI/WIOD_Projected/EnergyProjection_Aggregate.pdf",width=6,height=5)
ggplot() + 
  geom_line(data=df.WIOD[df.WIOD$Code=="TOTII_new",],aes(Year,`Total Energy Use (PJ)`,linetype="Data"),linewidth=1) +
  geom_line(data=df.WIOD.proj[df.WIOD.proj$Code=="TOTII_new",],aes(Year,`Total Energy Use (PJ)`,linetype="Projection"),linewidth=1) +
  theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + xlab("Year") + ylab("Total Energy Use (PJ)") + 
  ggtitle(paste0("Energy Use - ","Aggregate")) +
  scale_linetype_manual(name="Source",values=c("Data"="solid","Projection"="dashed")) 
dev.off()

pdf("TablesAndFigures/SI/WIOD_Projected/EmissionsProjection_Aggregate.pdf",width=6,height=5)
ggplot() + 
  geom_line(data=df.WIOD[df.WIOD$Code=="TOTII_new",],aes(Year,`CO2 Emissions (kt)`,linetype="Data"),linewidth=1) +
  geom_line(data=df.WIOD.proj[df.WIOD.proj$Code=="TOTII_new",],aes(Year,`CO2 Emissions (kt)`,linetype="Projection"),linewidth=1) +
  theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + xlab("Year") + ylab("CO2 Emissions (kt)") + 
  ggtitle(paste0("Emissions - ","Aggregate")) +
  scale_linetype_manual(name="Source",values=c("Data"="solid","Projection"="dashed")) 
dev.off()

pdf("TablesAndFigures/SI/WIOD_Projected/EnergyIntensityProjection_Aggregate.pdf",width=6,height=5)
ggplot() + 
  geom_line(data=df.WIOD[df.WIOD$Code=="TOTII_new",],aes(Year,`Energy Intensity (PJ/$BB)`,linetype="Data"),linewidth=1) +
  geom_line(data=df.WIOD.proj[df.WIOD.proj$Code=="TOTII_new",],aes(Year,`Energy Intensity (PJ/$BB)`,linetype="Projection"),linewidth=1) +
  theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + xlab("Year") + ylab("Energy Intensity (PJ/$BB)") + 
  ggtitle(paste0("Energy Intensity - ","Aggregate")) +
  scale_linetype_manual(name="Source",values=c("Data"="solid","Projection"="dashed")) 
dev.off()


pdf("TablesAndFigures/SI/WIOD_Projected/EmissionsIntensityProjection_Aggregate.pdf",width=6,height=5)
ggplot() + 
  geom_line(data=df.WIOD[df.WIOD$Code=="TOTII_new",],aes(Year,`Emissions Intensity (kt/PJ)`,linetype="Data"),linewidth=1) +
  geom_line(data=df.WIOD.proj[df.WIOD.proj$Code=="TOTII_new",],aes(Year,`Emissions Intensity (kt/PJ)`,linetype="Projection"),linewidth=1) +
  theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + xlab("Year") + ylab(expression("Emissions Intensity (ktC"*O[2]*"/PJ)")) + 
  ggtitle(paste0("Emissions Intensity - ","Aggregate")) +
  scale_linetype_manual(name="Source",values=c("Data"="solid","Projection"="dashed")) 
dev.off()

for(code in unique(df.WIOD[df.WIOD$Code!="U" & df.WIOD$Code!="TOTII_new","Code"])){
  
  
  Oplot <- ggplot() + 
    geom_line(data=df.WIOD[df.WIOD$Code==code,],aes(Year,`Output ($BB)`,linetype="Data"),linewidth=1) +
    geom_line(data=df.WIOD.proj[df.WIOD.proj$Code==code,],aes(Year,`Output ($BB)`,linetype="Projection"),linewidth=1) +
    theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + xlab("Year") + ylab("Output ($BB)") + 
    ggtitle(paste0("Output - ",code)) +
    scale_linetype_manual(name="Source",values=c("Data"="solid","Projection"="dashed")) 
  
  pdf(paste0("TablesAndFigures/SI/WIOD_Projected/OutputProjection_",code,".pdf"),width=6,height=5)
  print(Oplot)
  dev.off()
  
  
  Enplot <- ggplot() + 
    geom_line(data=df.WIOD[df.WIOD$Code==code,],aes(Year,`Total Energy Use (PJ)`,linetype="Data"),linewidth=1) +
    geom_line(data=df.WIOD.proj[df.WIOD.proj$Code==code,],aes(Year,`Total Energy Use (PJ)`,linetype="Projection"),linewidth=1) +
    theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + xlab("Year") + ylab("Total Energy Use (PJ)") + 
    ggtitle(paste0("Energy Use - ",code)) +
    scale_linetype_manual(name="Source",values=c("Data"="solid","Projection"="dashed")) 
  
  pdf(paste0("TablesAndFigures/SI/WIOD_Projected/EnergyProjection_",code,".pdf"),width=6,height=5)
  print(Enplot)
  dev.off()
  
  
  Emplot <- ggplot() + 
    geom_line(data=df.WIOD[df.WIOD$Code==code,],aes(Year,`CO2 Emissions (kt)`,linetype="Data"),linewidth=1) +
    geom_line(data=df.WIOD.proj[df.WIOD.proj$Code==code,],aes(Year,`CO2 Emissions (kt)`,linetype="Projection"),linewidth=1) +
    theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + xlab("Year") + ylab("CO2 Emissions (kt)") + 
    ggtitle(paste0("Emissions - ",code)) +
    scale_linetype_manual(name="Source",values=c("Data"="solid","Projection"="dashed")) 
  
  pdf(paste0("TablesAndFigures/SI/WIOD_Projected/EmissionsProjection_",code,".pdf"),width=6,height=5)
  print(Emplot)
  dev.off()
  
  
  
  EIplot <- ggplot() + 
    geom_line(data=df.WIOD[df.WIOD$Code==code,],aes(Year,`Energy Intensity (PJ/$BB)`,linetype="Data"),linewidth=1) +
    geom_line(data=df.WIOD.proj[df.WIOD.proj$Code==code,],aes(Year,`Energy Intensity (PJ/$BB)`,linetype="Projection"),linewidth=1) +
    theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + xlab("Year") + ylab("Energy Intensity (PJ/$BB)") + 
    ggtitle(paste0("Energy Intensity - ",code)) +
    scale_linetype_manual(name="Source",values=c("Data"="solid","Projection"="dashed")) 
  
  pdf(paste0("TablesAndFigures/SI/WIOD_Projected/EnergyIntensityProjection_",code,".pdf"),width=6,height=5)
  print(EIplot)
  dev.off()
  
  
  CIplot<- ggplot() + 
    geom_line(data=df.WIOD[df.WIOD$Code==code,],aes(Year,`Emissions Intensity (kt/PJ)`,linetype="Data"),linewidth=1) +
    geom_line(data=df.WIOD.proj[df.WIOD.proj$Code==code,],aes(Year,`Emissions Intensity (kt/PJ)`,linetype="Projection"),linewidth=1) +
    theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + xlab("Year") + ylab(expression("Emissions Intensity (ktC"*O[2]*"/PJ)")) + 
    ggtitle(paste0("Emissions Intensity - ",code)) +
    scale_linetype_manual(name="Source",values=c("Data"="solid","Projection"="dashed")) 
  
  pdf(paste0("TablesAndFigures/SI/WIOD_Projected/EmissionsIntensityProjection_",code,".pdf"),width=6,height=5)
  print(CIplot)
  dev.off()
  
}




################################################################################
################################################################################
# 3. Merge NAICS exposure weights (4-digit) with WIOD ISIC classifications
################################################################################
################################################################################

#############
# Create concordance array between ISIC codes and NAICS
#############

# Load NAICS 4-digit exposure weight data
df.NAICS.exp <- read.csv(file="Data/Output/NAICS4_Occupation_wagebillweighted_exposure.csv")

# Initialize array to store mapping
# Get 2-digit ISIC codes
df.ISIC_codes <- read_excel("Data/Input/CL_ACTIVITY_ISIC4_1.0.xls")
ISIC_codes_2_withletter <- array(unlist(df.ISIC_codes[df.ISIC_codes$`Annotation 1 - Hierarchical Level`=="Hierarchical level 2","Code Value"]))
ISIC_codes_2 <- substr(ISIC_codes_2_withletter,2,3)
# Get 4-digit NAICS codes
NAICS_codes_4 <- substr(unique(df.NAICS.exp$NAICS),1,4)
# Some NAICS codes are either multiples of 4 digit codes (e.g. ending in A1) or have bad match. For some, aggregate to 3 digit when possible. For others, aggregate to 2 digit. For others, no match.
missmatch3 <- c("325","327","332","333","337","423","424","445","484","522","523","531","532")
missmatch2 <- c("45","51")
# Create array using all matches
concord.arr.all <- array(0,dim=c(length(ISIC_codes_2),length(NAICS_codes_4)),dimnames=list(ISIC_codes_2,NAICS_codes_4))

# Fill in array using all matches
missing.naics = c()
for(naics in NAICS_codes_4){
  naics_lookup = naics
  
  # Replace with 2 or 3 digit naics when needed
  if(substr(naics,1,3)%in%missmatch3){naics_lookup = substr(naics,1,3)}
  if(substr(naics,1,2)%in%missmatch2){naics_lookup = substr(naics,1,2)}
  
  # Retreive isic code from concordance
  isic = concord_naics_isic(naics_lookup,"NAICS2017","ISIC4",dest.digit=2,all = TRUE)
  
  if(is.na(isic[[1]]$match[1])){
    print(paste0("No match found for ",naics))
    missing.naics = c(missing.naics,naics)
  }else{
    for(m in c(1:length(isic[[1]]$match))){
      concord.arr.all[isic[[1]]$match[m],naics] = concord.arr.all[isic[[1]]$match[m],naics] + isic[[1]]$weight[m]
    }
  }
}

# Create array using best match
concord.arr.best <- array(0,dim=c(length(ISIC_codes_2),length(NAICS_codes_4)),dimnames=list(ISIC_codes_2,NAICS_codes_4))

# Fill in array using best matches
missing.naics = c()
for(naics in NAICS_codes_4){
  naics_lookup = naics
  
  # Replace with 2 or 3 digit naics when needed
  if(substr(naics,1,3)%in%missmatch3){naics_lookup = substr(naics,1,3)}
  if(substr(naics,1,2)%in%missmatch2){naics_lookup = substr(naics,1,2)}
  
  # Retreive isic code from concordance
  isic = concord_naics_isic(naics_lookup,"NAICS2017","ISIC4",dest.digit=2,all = FALSE)
  
  if(is.na(isic)){ # If no match at 2-digit, throw warning
    missing.naics = c(missing.naics,naics)
    print(paste0("No match found for ",naics))
  }else{
    concord.arr.best[isic,naics] = 1
  }
}


#############
# Aggregate concordance array along ISIC codes to match WIOD
#############


# Data.frame matching ISIC codes in WIOD to all 2-digit ISIC codes
wiod_to_isic <- data.frame("ISIC2" = ISIC_codes_2,"WIOD_code" = NA)
letter_look_up = function(letter){
  nums = substr(ISIC_codes_2_withletter[substr(ISIC_codes_2_withletter,1,1)==letter],2,3)
  return(nums)
}
for(i in unique(df.WIOD$Code)[1:56]){
  if(nchar(i)==1){ # If WIOD code is a number
    nums = letter_look_up(i)
    for(n in nums){
      wiod_to_isic[wiod_to_isic$ISIC2==n,"WIOD_code"] = i
    }
  }else if(grepl("_",i)){# If WIOD code spans multiple codes
    if(i=="R_S"){# If R/S
      i_subs = c("R","S")
      for(is in i_subs){
        nums = letter_look_up(is)
        for(n in nums){
          wiod_to_isic[wiod_to_isic$ISIC2==n,"WIOD_code"] = i
        }
      }
    }else{# If multiple numbers
      nums = as.character(c(as.numeric(substr(i,2,3)):as.numeric(substr(i,5,6))))
      for(n in nums){
        wiod_to_isic[wiod_to_isic$ISIC2==n,"WIOD_code"] = i
      }
    }
  }else if(nchar(i)==3){
    wiod_to_isic[wiod_to_isic$ISIC2==substr(i,2,3),"WIOD_code"] = i
  }
}

# Aggregate along WIOD codes
concord.arr.all.WIOD <- concord.arr.best.WIOD <- array(0,dim=c(length(unique(df.WIOD$Code)[1:56]),length(NAICS_codes_4)),dimnames=list(unique(df.WIOD$Code)[1:56],NAICS_codes_4))
for(i in dimnames(concord.arr.best.WIOD)[[1]]){
  isics = wiod_to_isic[wiod_to_isic$WIOD_code==i,"ISIC2"]
  if(length(isics)==1){ # If only one 2-digit code in WIOD ISIC code, look up code
    concord.arr.all.WIOD[i,] <- concord.arr.all[isics,]
    concord.arr.best.WIOD[i,] <- concord.arr.best[isics,]
  }else{ # If multiple 2-digit code in WIOD ISIC code, aggregate over ISIC codes
    concord.arr.all.WIOD[i,] <- apply(concord.arr.all[isics,],2,sum)
    concord.arr.best.WIOD[i,] <- apply(concord.arr.best[isics,],2,sum)
  }
}




##########
# Merge NAICS data with WIOD
##########
# Get list of exposure metrics in NAICS dataset
expmeasures = colnames(df.NAICS.exp)[substr(colnames(df.NAICS.exp),1,8)=="expshare"]

# Duplicate WIOD data for merging
df.WIOD.exp.all <- df.WIOD.exp.best <- df.WIOD.proj
df.WIOD.exp.all[expmeasures] <- df.WIOD.exp.best[expmeasures] <- 0
for(naic in unique(df.NAICS.exp$NAICS)){
  if(sum(concord.arr.best.WIOD[,substr(naic,1,4)])==0){next} # If NAICS doesn't match, skip
  
  # Fill in for all match
  cells.all = names(which(concord.arr.all.WIOD[,substr(naic,1,4)]>0))
  weights.all = concord.arr.all.WIOD[cells.all,substr(naic,1,4)]
  if(length(cells.all)>1){ # If NAICS matches to multiple isic codes, weight using weights from concord function
    for(j in c(1:length(cells.all))){
      for(i in expmeasures){
        df.WIOD.exp.all[df.WIOD.exp.all$Code==cells.all[j],i] = df.WIOD.exp.all[df.WIOD.exp.all$Code==cells.all[j],i] + df.NAICS.exp[df.NAICS.exp$NAICS==naic,i]*weights.all[j]
      }
    }
  }else if(length(cells.all)==1){
    for(i in expmeasures){
      df.WIOD.exp.all[df.WIOD.exp.all$Code==cells.all,i] = df.WIOD.exp.all[df.WIOD.exp.all$Code==cells.all,i] + df.NAICS.exp[df.NAICS.exp$NAICS==naic,i]
    }
  }
  
  # Fill in for best match
  cells.best = names(which(concord.arr.best.WIOD[,substr(naic,1,4)]>0))
  weights.best = concord.arr.best.WIOD[cells.best,substr(naic,1,4)]
  for(i in expmeasures){
    df.WIOD.exp.best[df.WIOD.exp.best$Code==cells.best,i] = df.WIOD.exp.best[df.WIOD.exp.best$Code==cells.best,i] + df.NAICS.exp[df.NAICS.exp$NAICS==naic,i]
  }
}







################################################################################
################################################################################
# 4. Make Exposure plots
################################################################################
################################################################################

# Set plot data (Projected)
df.WIOD.plot = df.WIOD.exp.best[df.WIOD.exp.best$Year %in% c(2015:2023) & df.WIOD.exp.best$Code!="U" & df.WIOD.exp.best$Code!="TOTII_new",]

# Average over projected years
df.WIOD.plot <- df.WIOD.plot %>%
  group_by(Code,Description) %>%
  summarize_all(mean,na.rm=T)



# Adjust to be share of VA_k rather than relative to GDP
for(i in expmeasures[c(2,4,6,7,10,12,14,16)]){
  df.WIOD.plot[i] = df.WIOD.plot[i]*23e9/(df.WIOD.plot['Output ($BB)']*1e6)*0.23*0.27/10
}

df.WIOD.plot = df.WIOD.plot[order(df.WIOD.plot$expshare_automation_binary_laborsharegdp_weight),]
df.WIOD.plot['Code_f'] = factor(df.WIOD.plot$Code,levels=df.WIOD.plot$Code)
df.WIOD.plot['Description_short'] = df.WIOD.plot['Description']
df.WIOD.plot['Description_f2'] = 0
for(i in c(1:dim(df.WIOD.plot)[1])){
  df.WIOD.plot[i,'Description_f2'] = i
  if(nchar(df.WIOD.plot[i,"Description_short"])>43) {
    df.WIOD.plot[i,"Description_short"] = paste0(substr(df.WIOD.plot[i,"Description_short"],1,40),"...")
  }
  
}
df.WIOD.plot['Description_f'] = factor(df.WIOD.plot$Description_short,levels=df.WIOD.plot$Description_short)




pdf("TablesAndFigures/SI/WIOD_Projected/ISICRev4_Industry_Wagebillweighted_Exposure.pdf",width=3,height=5)
ggplot(df.WIOD.plot[order(df.WIOD.plot$expshare_automation_binary_laborsharegdp_weight),]) +
  geom_bar(aes(x=expshare_automation_binary_laborsharegdp_weight*100,y=Description_f,fill=Description_f2),stat='identity',width=.5) +
  theme_classic() + theme(axis.text.y = element_text(size = 5),axis.text.x = element_text(size = 5),text = element_text(size=8)) +
  scale_fill_gradient(low = "lightgreen", high = "darkgreen") +
  theme(legend.position = "none") +
  scale_x_continuous(position = "top") +
  xlab(expression(frac(Delta*A[k],A[k])*" (%)")) + ylab("ISIC Rev.4 Industry") #+ ggtitle("Industry wage-bill weighted exposure shares (ISIC Rev.4)")
dev.off()



################################################################################
################################################################################
# 5. Make Energy and Emissions change plots (Aggregate)
################################################################################
################################################################################
# Set exposure score to be used
expscore = "automation_binary"
expscore.weight = "_laborsharegdp_weight"


df.WIOD.plot['Change in GDP'] = (df.WIOD.plot[paste0('expshare_',expscore,expscore.weight)])*df.WIOD.plot$`Output ($BB)`
df.WIOD.plot['Change in Energy'] = df.WIOD.plot['Change in GDP']*df.WIOD.plot$`Energy Intensity (PJ/$BB)`
df.WIOD.plot['Change in Emissions'] = df.WIOD.plot['Change in Energy']*df.WIOD.plot['Emissions Intensity (kt/PJ)']
df.WIOD.plot['Change in GDP (upper)'] = (df.WIOD.plot[paste0('expshare_',expscore,"_upper",expscore.weight)])*df.WIOD.plot$`Output ($BB)`
df.WIOD.plot['Change in Energy (upper)'] = df.WIOD.plot['Change in GDP (upper)']*df.WIOD.plot$`Energy Intensity (PJ/$BB)`
df.WIOD.plot['Change in Emissions (upper)'] = df.WIOD.plot['Change in Energy (upper)']*df.WIOD.plot['Emissions Intensity (kt/PJ)']
df.WIOD.plot['Change in GDP (lower)'] = (df.WIOD.plot[paste0('expshare_',expscore,"_lower",expscore.weight)])*df.WIOD.plot$`Output ($BB)`
df.WIOD.plot['Change in Energy (lower)'] = df.WIOD.plot['Change in GDP (lower)']*df.WIOD.plot$`Energy Intensity (PJ/$BB)`
df.WIOD.plot['Change in Emissions (lower)'] = df.WIOD.plot['Change in Energy (lower)']*df.WIOD.plot['Emissions Intensity (kt/PJ)']


print(paste0("Labor share GDP weighted exposure (aggregate) = ",as.character(sum((df.WIOD.plot[paste0('expshare_',expscore,"_upper",expscore.weight)])))))
# print(paste0("Change in energy use (aggregate) = ",as.character(round(sum((df.WIOD.plot[paste0('expshare_',expscore,expscore.weight)]*0.23*0.27)/10)*df.WIOD[df.WIOD$Code=="TOTII_new" & df.WIOD$Year==2014,"Output ($BB)"]*df.WIOD[df.WIOD$Code=="TOTII_new" & df.WIOD$Year==2014,"Energy Intensity (PJ/$BB)"]),3)))
print(paste0("Change in GDP (by industry) = ",as.character(round(sum(df.WIOD.plot$`Change in GDP`),3))))
print(paste0("Change in energy use (by industry) = ",as.character(round(sum(df.WIOD.plot$`Change in Energy`),3))))
print(paste0("Change in emissions (by industry) = ",as.character(round(sum(df.WIOD.plot$`Change in Emissions`),3))))

# Export relevant data as table
write.csv(df.WIOD.plot,"TablesAndFigures/SI/WIOD_Projected/FinalEstimates.csv",row.names=F)

# Set codes to highlight (Publishing activities, Education,Wholesale trade and repair of motor vehicles)
codesofinterest = c("J58","P85","G45")
codesofinterest.color <- c("J58" = "red", "P85" = "blue", "G45" = "green")



# Create a data frame

xvar = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Output ($BB)"]
yvar_min = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"expshare_automation_binary_lower_laborsharegdp_weight"]
yvar_max = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"expshare_automation_binary_upper_laborsharegdp_weight"]
x_range = seq(min(xvar),max(xvar),(max(xvar)-min(xvar))/100)
df001 = data.frame("x"=x_range,"y"=1e-5/(x_range))
df10 = data.frame("x"=x_range,"y"=1e-3/(x_range))
df100 = data.frame("x"=x_range,"y"=1e-1/(x_range))
df10000 = data.frame("x"=x_range,"y"=1e1/(x_range))
gdp = sum(df.WIOD.plot$`Output ($BB)`)
plt.ratio = 2/5
x.pos = exp(log(min(x_range)+0.00025*(max(x_range)-min(x_range))))

plt1 <- ggplot() +
  geom_point(data=df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,],aes(`Output ($BB)`,expshare_automation_binary_laborsharegdp_weight),color="grey") +
  geom_point(data=df.WIOD.plot[df.WIOD.plot$Code %in% codesofinterest,],aes(`Output ($BB)`,expshare_automation_binary_laborsharegdp_weight,color=Code)) +
  geom_linerange(data=df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,],aes(x=`Output ($BB)`,ymin=expshare_automation_binary_lower_laborsharegdp_weight,ymax=expshare_automation_binary_upper_laborsharegdp_weight),color="grey") +
  geom_linerange(data=df.WIOD.plot[df.WIOD.plot$Code %in% codesofinterest,],aes(x=`Output ($BB)`,ymin=expshare_automation_binary_lower_laborsharegdp_weight,ymax=expshare_automation_binary_upper_laborsharegdp_weight,color=Code)) +
  scale_x_continuous(trans='log10') + scale_y_continuous(trans='log10',limits=c(1e-7,2e-2)) +
  geom_line(data = df001,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=1e-6,label="$10 Thousand",angle=180/pi*atan(-1*plt.ratio),hjust = 0) +
  geom_line(data = df10,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=1e-4,label="$1 Million",angle=180/pi*atan(-1*plt.ratio),hjust = 0) +
  geom_line(data = df100,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=1e-2,label="$100 Million",angle=180/pi*atan(-1*plt.ratio),hjust = 0) +
  geom_line(data = df10000,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=750,y=2e-2,label="$10 Billion",angle=180/pi*atan(-1*plt.ratio),hjust = 0) +
  scale_color_manual(values = codesofinterest.color) +
  ylab(expression(Delta*A[k]*"/"*A[k])) + xlab("Industry Size ($BB)") +
  theme_classic() + theme(legend.position="none") #+ coord_fixed(ratio=plt.ratio)

print(plt1)


pdf("TablesAndFigures/SI/WIOD_Projected/ShockvsValueAdded.pdf",width=5,height=5)
print(plt1)
dev.off()


xvar = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Energy Intensity (PJ/$BB)"]
yvar_min = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Change in GDP (lower)"]
yvar_max = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Change in GDP (upper)"]
x_range = seq(min(xvar),max(xvar),(max(xvar)-min(xvar))/1000)
x_range = x_range[x_range>0]
df0001 = data.frame("x"=x_range,"y"=1e-5/x_range)
df10 = data.frame("x"=x_range,"y"=1e-2/x_range)
df100 = data.frame("x"=x_range,"y"=1/x_range)
df10000 = data.frame("x"=x_range,"y"=1e2/x_range)
plt.ratio = 2/3
x.pos = exp(log(min(x_range)+0.0000025*(max(x_range)-min(x_range))))

plt2 <- ggplot() +
  geom_point(data=df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,],aes((`Energy Intensity (PJ/$BB)`),`Change in GDP`),color="grey") +
  geom_point(data=df.WIOD.plot[df.WIOD.plot$Code %in% codesofinterest,],aes((`Energy Intensity (PJ/$BB)`),`Change in GDP`,color=Code)) +
  geom_linerange(data=df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,],aes(x=(`Energy Intensity (PJ/$BB)`),ymin=`Change in GDP (lower)`,ymax=`Change in GDP (upper)`),color="grey") +
  geom_linerange(data=df.WIOD.plot[df.WIOD.plot$Code %in% codesofinterest,],aes(x=(`Energy Intensity (PJ/$BB)`),ymin=`Change in GDP (lower)`,ymax=`Change in GDP (upper)`,color=Code)) +
  geom_line(data = df0001,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=3e-3,label="10 GJ",angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  geom_line(data = df10,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=3,label="10 TJ",angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  geom_line(data = df100,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=1.1e-1,y=1.5e1,label="1 PJ",angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  geom_line(data = df10000,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=0.8e1,y=2e1,label="100 PJ",angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  scale_color_manual(values = codesofinterest.color) +
  scale_x_continuous(trans='log10') + scale_y_continuous(trans='log10',breaks=c(1e-5,1e-3,1e-1,1e1),limits=c(0.5e-5,2e1)) +#,breaks=c(1e-3,1,1e3),limits = c(1e-3,1e6)) +
  ylab("Change in GDP ($BB)") + xlab("Energy Intensity (PJ/$BB)") +
  theme_classic() + theme(legend.position="none") #+ coord_fixed(ratio=plt.ratio)

print(plt2)

pdf("TablesAndFigures/SI/WIOD_Projected/ChangeinGDPvsEnergyIntensity.pdf",width=5,height=5)
print(plt2)
dev.off()

xvar = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Emissions Intensity (kt/PJ)"]
yvar_min = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Change in Energy (lower)"]
yvar_max = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Change in Energy (upper)"]
x_range = seq(min(xvar),max(xvar),(max(xvar)-min(xvar))/100)
y_range = 1/x_range
df00001 = data.frame("x"=x_range,"y"=1e-3/x_range)
df1 = data.frame("x"=x_range,"y"=1e-1/x_range)
df100 = data.frame("x"=x_range,"y"=1e1/x_range)
df1000 = data.frame("x"=x_range,"y"=1e3/x_range)
plt.ratio = 2/8
x.pos = exp(log(min(x_range)+0.0000025*(max(x_range)-min(x_range))))

plt3 <- ggplot() +
  geom_point(data=df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,],aes(`Emissions Intensity (kt/PJ)`,`Change in Energy`),color="grey") +
  geom_point(data=df.WIOD.plot[df.WIOD.plot$Code %in% codesofinterest,],aes(`Emissions Intensity (kt/PJ)`,`Change in Energy`,color=Code)) +
  geom_linerange(data=df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,],aes(x=(`Emissions Intensity (kt/PJ)`),ymin=`Change in Energy (lower)`,ymax=`Change in Energy (upper)`),color="grey") +
  geom_linerange(data=df.WIOD.plot[df.WIOD.plot$Code %in% codesofinterest,],aes(x=(`Emissions Intensity (kt/PJ)`),ymin=`Change in Energy (lower)`,ymax=`Change in Energy (upper)`,color=Code)) +
  scale_x_continuous(trans='log10') + scale_y_continuous(trans='log10',limits = c(0.5e-5,2e2),breaks=c(1e-5,1e-2,1,1e2)) + #,breaks=c(1e-3,1,1e3),limits = c(1e-5,1e5)) +
  geom_line(data = df00001,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=0.5e-3,label=expression("1 tC"*O[2]),angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  geom_line(data = df1,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=0.5e-1,label=expression("100 tC"*O[2]),angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  geom_line(data = df100,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=0.5e1,label=expression("10 ktC"*O[2]),angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  geom_line(data = df1000,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=9,y=2e2,label=expression("1 MtC"*O[2]),angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  scale_color_manual(values = codesofinterest.color) +
  ylab("Change in Energy (PJ)") + xlab("Emission Intensity (ktCO2/PJ)") +
  theme_classic() + theme(legend.position="none") #+ coord_fixed(ratio=plt.ratio)

print(plt3)


pdf("TablesAndFigures/SI/WIOD_Projected/ChangeinEnergyvsEmissionsIntensity.pdf",width=5,height=5)
print(plt3)
dev.off()


grid.arrange(plt1, plt2, plt3, ncol=3)





################################################################################
################################################################################
# 6. Make Energy and Emissions change plots (By Industry)
################################################################################
################################################################################
df.WIOD.plot <- df.WIOD.plot[order(df.WIOD.plot$Code),]
df.WIOD.plot['Description_f'] = factor(df.WIOD.plot$Description_short,levels=df.WIOD.plot$Description_short)

# Plot impact on GDP by Industry
pdf("TablesAndFigures/SI/WIOD_Projected/Industry_GDP_Impact.pdf",width=3,height=5)
ggplot() +
  geom_bar(data=df.WIOD.plot[order(df.WIOD.plot$Code),],aes(x=`Change in GDP`,y=Description_f,fill=Description_f2),stat='identity',width=.5) +
  theme_classic() + theme(axis.text.y = element_text(size = 5),axis.text.x = element_text(size = 5),text = element_text(size=8)) +
  scale_fill_gradient(low = "lightgreen", high = "darkgreen") +
  theme(legend.position = "none") +
  scale_x_continuous(position = "top") +
  xlab("Change in GDP ($BB)") + ylab("ISIC Rev.4 Industry") 
dev.off()


# Plot impact on Energy Use by Industry
pdf("TablesAndFigures/SI/WIOD_Projected/Industry_EnergyUse_Impact.pdf",width=3,height=5)
ggplot() +
  geom_bar(data=df.WIOD.plot,aes(x=`Change in Energy`,y=Description_f,fill=Description_f2),stat='identity',width=.5) +
  theme_classic() + theme(axis.text.y = element_text(size = 5),axis.text.x = element_text(size = 5),text = element_text(size=8)) +
  scale_fill_gradient(low = "lightgreen", high = "darkgreen") +
  theme(legend.position = "none") +
  scale_x_continuous(position = "top") +
  xlab("Change in Energy (PJ)") + ylab("ISIC Rev.4 Industry") #+ ggtitle("Industry wage-bill weighted exposure shares (ISIC Rev.4)")
dev.off()


# Plot impact on Emissions by Industry
pdf("TablesAndFigures/SI/WIOD_Projected/Industry_Emissions_Impact.pdf",width=3,height=5)
ggplot() +
  geom_bar(data=df.WIOD.plot,aes(x=`Change in Emissions`,y=Description_f,fill=Description_f2),stat='identity',width=.5) +
  theme_classic() + theme(axis.text.y = element_text(size = 5),axis.text.x = element_text(size = 5),text = element_text(size=8)) +
  scale_fill_gradient(low = "lightgreen", high = "darkgreen") +
  theme(legend.position = "none") +
  scale_x_continuous(position = "top") +
  xlab("Change in Emissions (ktCO2)") + ylab("ISIC Rev.4 Industry") #+ ggtitle("Industry wage-bill weighted exposure shares (ISIC Rev.4)")
dev.off()









